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Abstract. We study the off-equilibrium relaxational dynamics at criticality in the 
three-dimensional Blume-Capel model whose static critical behaviour belongs to the 
3d-Ising universality class. Using "improved" Hamiltonian (the leading corrections 
to scaling have vanishing amplitude) we perform Monte Carlo simulations of the 
relaxational dynamics after a quench from T = oo to Tc- Analysing the off-equilibrium 
dynamics at Tc we obtain an estimate of the dynamical critical exponent z = 2.020(8) 
that is perfectly consistent with the Field Theory predictions. 



1. Introduction 

Many statistical and dynamical properties of different physical systems show 
"anomalies" if the thermodynamics parameters approach a critical point. This singular 
behaviour exhibits a great universality. In fact, many properties of a system close 
to a continuos phase transition turn out to be largely independent of the microscopic 
details of the interactions between the individual atoms or molecules. They fall into 
one universality class characterized only by global features such as the symmetries 
of the underlying Hamiltonian, the number of spatial dimensions of the system, the 
number of components of order parameter and so on. Typically, close to a critical 
point, the correlation length and the other thermodynamics quantities exhibit power- 
law dependences on the parameters specifying the distance away from the critical point. 
These powers, or critical exponents, are pure numbers, usually not integers or even 
simple rational numbers, which depend only on the universality class [H El |3l H] . 

In critical dynamics we are interested in the universal features of the relaxational 
process: after preparing the system in a macroscopical initial state, we leave it to 
relax toward equilibrium. Upon approaching a critical point, the typical time scale 
of dynamics of the fluctuations around the equilibrium state diverges as ~ C,^ {critical 
slowing down), where ^ is the correlation length and z is the dynamical critical exponent. 

From theoretical point of view, after a critical quench, the dynamical behaviour far 
from equilibrium shows some universal features associated with the universal properties 
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of the critical point. Thus, studying the off-equihbrium critical dynamics it is possible 
to acquire informations about the static and the dynamic universality class of the model 
and to get estimates of many different critical quantities. 

So far, despite the efforts for understanding the off-equilibrium critical properties 
in many different systems [5], like two-dimensional Ising and Potts models [B], classical 
XY systems [3 |8], disordered quantum spin chains [9], randomly diluted classical 
models [TU], there is considerable theoretical and experimental interest to study some 
more fundamental systems, like 3d-Ising model, specially for having better numerical 
estimates of the critical exponents. 

Indeed, some differences in the determination of the dynamical critical exponent 
z for the three-dimensional Ising universality class yet exist between the numerical 
simulations, like Monte Carlo simulation at equilibrium [TTl [12] {z = 2.04(3), z = 
2.10(2)), damage spreading [El [H US] (z = 2.034(4), z = 2.04(1), z = 2.05(5)), critical 
relaxation [161 [H] (z = 2.05(2), z = 2.04(2)), short-time dynamics [H] (z = 2.042(6)), 
and the Field Theory predictions [IHl [20] {z ~ 2.02, z ~ 2.017). The non-equilibrium 
relaxation method was also used to test the dynamic universality hypothesis in three- 
dimensional Ising models and two-dimensional three-states Potts models [21] . 

In order to have a better estimate of the dynamical critical exponent, we have 
studied a purely relaxational dynamics without conserved order parameter (also known 
as model A[|]) in 3d-Ising system; in particular, we have focused on the growth of 
correlations during the first instants of the evolution of the system after a quench 
from T = oo to T = Tc (off-equilibrium critical dynamics). We have used a Blume- 
Capel Hamiltonian with improved parameters (see Sec. 2 for definitions and details). 
According to static and dynamics TZQ (Renormalization Group), the critical behaviour 
of such a system belongs into the Ising model universality class. For that reason, the 
dynamical critical exponent z is the same. In addition, the leading scaling corrections 
appearing in dynamical quantities are governed by the same TZQ operators that control 
the non-asymptotic behaviour of static quantities and thus, are characterized by the 
same exponents as in the static case, i.e., by Ui = 0.84(4) and W2 = 1.67(11) [11I221I23]- 
As a consequence, the equilibrium improved Hamiltonians will cancel the leading scaling 
corrections also in the dynamical quantities. Therefore, the most precise estimates of 
dynamical universal quantities should be obtained in improved models, as in the static 
case [2i[[25]. 

I The model A dynamical universality class can be specified in terms of the stochastic Langevin 
equation 

5t<^(x,t) = -r!/^+C(a;,i), 
oip(x, t) 

where ip{x,t) is the order parameter field, 'H[</?] is the static reduced Hamiltonian, f2 is a kinetic 
coefficient and Ci^it) is a zero-mean Gaussian white noise with correlations 

{ax,my,s)) =2nS{x-y)S{t-s). 

The critical dynamics of some anisotropic magnets and alloys are described by the previous dynamical 
equations with H given by the effective Hamiltonian of the Ising universality class. 
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We have been also careful at new scaling corrections introduced by dynamics. 
We have seen that the out-of-equilibrium dynamical quantities depend in general on 
the initial conditions and show analytical corrections tuned by Tq"^ that controls the 
correlations of the order parameter at initial time surface j§] Using the dynamical TZQ 
methods one finds that it is possible to forget such as "dynamical" corrections in the 
data analysis. 

We have performed MC (Monte Carlo) simulations for different lattice size L and 
we have observed the initial temporal behaviour of susceptibility x{^) L) and correlation 
length C,{t,L). We have noticed that for L/C,{t,L) > 5 the finite-size effects are small 
and the behaviour of x(t, L) is independent from lattice length. Using the susceptibility 
behaviour in this range, we have got our best value for the dynamical critical exponent 

z = 2.020(8). (1) 

Our result is perfectly in agreement with FT (Field Theory) approach, both e-expansion 
[19] (z ^ 2.02) and fixed dimension [20] (z ^ 2.017). 

The paper is organized as follows. In Sec. 2 we define the models and the 
observables, giving also the scaling form predicted by the TZQ analysis. In Sec. 3 
we report the analysis of MC simulations of the off-equilibrium relaxational critical 
behaviour in a quench from T = cxd to T^. Finally, we draw our conclusions in Sec. 4. 
Some details on the MC algorithm are discussed in Appendix A. 

2. Model and Observables 

We consider the Blume-Capel model with Hamiltonian 

H = -pJ2''^'^ + ^Yl''^ (2) 

{id) i 

where the Sj variables take the values {— 1,0, +1}, and the i index runs over a cubic 
lattice with length L and periodic boundary conditions. The brackets {i,j) indicate 
that sum runs over neighbors sites. /3 and D parameters are tuned at their improved 
critical values [^ITTll^: 

D* = 0.641, (3c = 0.3856717. (3) 

According D parameter at the improved value D* minimizes the leading order 
scaling corrections and improves the dynamical behaviour also for small times. In 
short, using improved models we can obtain more reliable estimates about universal 
dynamical quantities. In fact, FSS (Finite Size Scaling) tells that, at criticality, in 
improved systems, the scaling corrections behave as L~^^ (with corresponding dynamic 
scaling corrections t~^'^^^), where U2 is the next to leading scaling corrections exponent; 
otherwise, in a generic system, the approach to the thermodynamics limit (L — )• oo) 
is slower, with scaling corrections that behave as L^'^^ (with corresponding dynamical 
behaviour t~'^^/^), with the leading scaling corrections exponent coi. 

§ If (p{x,t) is the order parameter, then tq is defined via ((^(a;, to)v(2/j ^o)) = '''o^^i^ ~ u)- 
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In this improved model we have measured the Susceptibihty 

i 

and the Dynamics Structur Factor C^it) at least impulse (|g| = l-a jV) 

C(t) = ^ + C^.(t) + (5) 

with 

C,it) = jr^\j2e''-'^s^\ (6) 

i 

where fi = {xi,yi,Zi) is the coordinate of the lattice site and i runs over all the cubic 
lattice. 

The system is invariant under axes permutation, so we have taken the average over 
the three spatial directions (x, y, z) in order to improve the statistical average. By ^ 
and ^ it is possible to define the Correlation Length ^ using the discrete form 



^^^^"V 4sin^(vr/L) " 
2.1. Susceptibility scaling form 

Using the dynamical TZQ methods and, taking into account further scaling corrections 
consequent by static properties of the Ising-like universality class, we expect that the 
magnetic susceptibility x{^)i the infinite- volume limit, grows with MC time t as [10] 

X{t) = x^tPF^{^i,t-^\^i^t-^\ . . . ,/i„t-^";/i.„r^^o), (8) 

where [27] 

p = T] = 0.0364(2). (9) 

Fy. is analytic function of arguments. Scaling fields /ij rule the statical scaling corrections, 
while iXt-q is another scaling field linked to the off-equilibrium dynamics. Taking into 
account only the first two statical scaling fields and expanding in Taylor series, we 
obtain: 

x{t) = xot^(l + Cnr"' + c^i2t-'"^ + • ■ • + C2it-^' + C22t-^''' + • • • 

+ Cj^^^t-^-o + o)t-2-^o +...). (10) 



According to the FT perturbative analysis, the leading scaling-corrections exponents 
should be the same as those that occur in equilibrium (static or dynamics) correlation 
functions. Therefore, we expect 

= y ~ 0.42(2), V2 = ^- 0.83(5), (11) 

where we have used = 0.84(4), uj2 = 1.67(11) « 2a;i and z = 2 [2 The 
leading scaling correction proportional to t~^^ (and also all corrections of the form t~^^^) 
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vanishes in improved models. Moreover, we discard the further scahng corrections due 
to Tq"^ parameter because they are lower than corrections introduced by V2 exponent, 
since = z and Vr^ = <^to/z = 1. 

Indeed, in the field-theoretical approach to non-equilibrium dynamics the 
expectation value of a generic observable has to be computed with a total dynamical 
functional given by S[(p, if] + T-Lo['^o], where 5* is the dynamical functional representation 
of the Langevin equation of the model A and 'Holv'o] = / d^XY'^lix) provides the initial 
conditions. Then, following the standard methods of the dynamical TZQ [29l |30j, the 
addition of 'Ho[v5o] breaks the TTI (Time-Translational Invariance) and gives rise to 
new divergences in perturbation theory whenever time approaches the 'time surface' 
located at to = 0. It is possible to remove this new singularities without introducing 
new renormalization constants. Finally, it is possible to show that the TZQ dimension of 
firo is exactly = z. For that reason, we set henceforth firo = 0. 

Equation (10) is valid only in the infinite-volume limit. For a finite system of size 
L we expect 

X{t, L) = xot'EoitL^ni^ + C2ir^'E,{tL-') + ■■■), (12) 

where we have set the improved conditions and forgotten Tq"^ corrections. Moreover, 
En{x) are universal functions satisfying 

^n(O) = 1, 

Eq{x) ~ x-P, En^x) ~ X ^ OO. (13) 

Even though, for improved models, from the theoretical point of view, it is perfectly 
correct to discard the leading-order scaling corrections, we should verify the absence of 
such a corrections in the data analysis. On the other hand, discarding corrections due 
to Tq"^ should not affect the data analysis. Indeed, for our data range, the difference 
between t~o-83 ^-i jg gj^all enough to allow us to absorb the effects of Tq"^ into the 
corrections introduced by V2 exponent (See Section 3 for more details). 



3. Monte Carlo Simulation 



The off-equilibrium relaxational dynamics of the Ising model has already been 
investigated in [18]. Their result does not agree with the above- reported predictions. 
Reference [IB] obtains z = 2.042(6). It is quite difficult to reconcile these results with the 
FT predictions; indeed, in [18] they do not take into account the corrections proportional 
to t~^''^'^, i.e. controlled by the leading scaling-corrections exponent Ui = 0.84(4), 
therefore that result could be affected by this assumption. 

Others numerical MC works [161 EZ] investigate the critical relaxation starting from 
ordered states. Although their estimates of the dynamical critical exponent are still 
compatible within the error bar with the theoretical predictions (obtaining z = 2.05(2) 
and z = 2.04(2)), we suppose that by using improved Hamiltonian we could obtain more 
reliable estimates of the dynamical critical exponent. 
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Figure 1. On the left, the effective exponent peff(t,L) for L = 48,64,96. On the 
right, the same exponent versus the rescaled time t' = i(_L/96)^^. 



Some effort was yet spent on the analysis of the 3d-Ising universahty class by using 
improved Hamiltonian |25] . Recently, the estimates of D* and /3c for the Blume- 
Capel Hamiltonian parameters was improved, providing also new estimates for critical 
exponents z/, rj and Ui |31j. Thus, in the following, we further investigate this issue. 
We study the Metropolis dynamics of the Ising universality class after a quench from 
T = oo to Tc- This represents a nontrivial check of the FT predictions. 



3.1. Estimate of the dynamic critical exponent z 

Since the model is improved, using only the next-to-leading scaling correction, we predict 



Xit) = Xot'il + C2it-^' + ■■■) 
We define an effective exponent 

ln[x{2t)/x{t)] 



Peffit) 



In 2 



which, for if: — )■ oo, behaves as 

Peffit) =P + Cpt-''^ + 



■"2 



In 2 



On a finite lattice, (16) is replaced by 



Peffit, L) =p + eo{tL 



(14) 



(15) 



(16) 



(17) 



where we have neglected large-t corrections and eo{x) is a universal function (apart from 
a normalization of the argument) such that 

eo(0) = 0, lim eo(x) 



-P- 



We have performed off-equilibrium MC simulations on cubic lattices with periodic 
boundary conditions and size L = 48,64,96 at /3 = 0.3856717 and D = 0.641 (our 
presently best estimate of D* is D = 0.641(8)). For each lattice we have averaged over 
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□ L = 48, t_max = 1 60 
o L = 64, t_max = 240 
_ A L = 64, t_max = 280 
« L = 96, t_max = 500 
V L = 96, t_max = 700 
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□ L = 48, t_max = 160 
o L = 64, t_max = 240 
A L = 64, t_max = 280 
o L = 96, t_max = 500 
V L = 96, t_max = 700 
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Figure 2. On the left, the resuhs of fits \nx{t) = Inxo + plnt + C2it""^ for 
imin ^ i ^ t^ax- The lines correspond to the estimate p = 0.972(3). On the right, the 
results of fits h\x{t) = Inxo + plni for data that do not show corrections to scaling. 
The lines correspond to the estimate p = 0.975(2). Some data are slightly shifted along 
the X axis to make them visible. 



2 ■ 10^ Markov chains. For each chain we have set different high-temperature initial 
conditions and, after a quench at critical point, we have done 400 lattice sweeps (700 
for L = 96) using the Metropolis algorithm described in Appendix A. 

In Figjl] (on the left) we show Peffit, L) for L = 48, 64, 96. It clearly shows finite- 
size corrections, and, for each L, peff{t,L) follows the infinite-volume curve up to an 
L-dependent value tmax{L). As shown by Figjl] (on the right), where Peffit, L) is plotted 
versus t^esc = ^(-^/96)~^, finite-size effects are consistent with (17). Thus, the value 
tmax{L), after which finite-size effects cannot be neglected, increases as L^. Infinite- 
volume quantities, such as peff{t), must be obtained from the data at t < tmax{L). Fig|l] 
indicates that, with the statistical errors of our data, tmax{L) ^ 90, 150 for L ^ 48, 64; 
moreover, using FSS Ansatz, we also predict tmax{L = 96) ~ 350|||] 

Since Peff{t, L) is defined using data at t and 2t, this implies that, for L = 96, only 
data corresponding to t < 700 have negligible finite-size effects within our error bars. 
Finite-size effects give rise to systematic error in the estimate of p. As is clear from 
FigjT| they yield smaller values of p, and therefore larger values of z. 

After selecting the temporal range without finite-size effects, from the susceptibility 
behaviour in this range, we have achieved our estimate of the dynamical critical exponent 
z. We have fitted the data by using (14): we have defined the temporal range [tmin, imax] 
with tmax{L = 48) = 160, t,nax{L = 64) = 240,280 and tmax{L = 96) = 500,700; for 
each tmax, tmin was moved in the interval [5, tmax/S]. In this way we have checked the 
stability of our result. 

In Figjl] (on the left) we plot the results of fits of In x(L, t) for tmin ^ t ^ tmax to 
Inxo + pint + C2it~^^ with V2 = 0.83, for L = 48,64,96. Finally, in Tabjl] we report 

II In Figjl] we have used our estimate for the dynamical critical exponent z. Instead, for the scaling 
analysis of tmax we have used the mean field approximation value z = 2. 
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our best estimate of p for different L and tmax getting into account the next-to-leading 
scaling corrections. Data corresponding to different lattice sizes are not correlated. 
Moreover, there are not finite-size effects in the data, thus we have got the estimate of 
the p exponent (in the thermodynamic limit L — )■ oo): 

p = 0.972(3), (19) 

and, using ([9]), we obtain 

z = 2.020(6). (20) 

Finally, in Figjs] we plot peffit.L) for L = 48,64,96 versus t'^^ ^j^h V2 = 0.83. 
Finite-size effects are negligible for t-""^ > tmaxiLy^ ^ 0.015,0.009,0.004, for L = 
48, 64, 96 respectively. The data satisfying this inequality clearly follow a unique curve, 
which is expected to behave as p + ct~'"'^ for sufficiently large values of t. The data 
plotted in Figjs] clearly show such a behaviour in the region t~'"^ < 0.2 (corresponding 
to t > 6), and are perfectly compatible with the value p = 0.972(3). 

Result without corrections to scaling - After that, we have also achieved an estimate 
of the exponent z by using data range that should not show next-to-leading scaling 
corrections. Thus, we have fitted the data with Inx(t) = Inxo + pint. The results are 
shown in Tab |2] and in Fig|2] (on the right). 

Also in this case, data corresponding to different lattice sizes are not correlated and 
there are not finite-size effects. We have got the estimate 

p = 0.975(2), z = 2.014(4), (21) 



perfectly consistent with (20). 



Errors analysis - As already mentioned at the end of Section 2, we can discard the 
corrections due to Tq'^ (~ t^^) since they are of the same order of those proportional to 
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Table 1. Some best estimate of p for different L, tmax and t, 



t We have seen, in fact, that fitting the data adding a new term proportional to 
did not change the results and introduced instability to the fit. It is easy to see that 
for t e [20,500] one has rO-^^ ^ [0.08,0.006] and G [0.05,0.002], thus f-^ ^ t'^-^^. 
As a consequence, the effects of are absorbed by 

On the other hand, for assuring that our data analysis was properly correct, we 
checked the previous results by adding, to the fit formula, the leading-order scaling- 
corrections term. In improved models such a term (~ t'""^) should be absent. We 
verified by fitting the data with the formula Inx(t) = Inxo + pint + Cut'""^ + 6*21^""^ 
and moving the fitting data range in the same intervals as we have explained in the 
previous paragraphs. Finally, we checked the behaviour of Cu (see figure |4]) and we 
verified his compatibility with zero value. In conclusion, taking into account the t~^^ 
term, we did not observe any appreciably changing in the estimate of the dynamical 
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Table 2. Some best estimate of p for different i, tmax and tmin- Here we only consider 
data that do not show corrections to scahng. 

critical exponent z. Furthermore, checking the fit by using the least-square test, we 
have not noticed any improvement by respect to the previous fits. 

Furthermore, like in others works that used improved models [I0l[32], we have took 
care of further corrections by checking the stability of the fits by varying the temporal 
range in which the fits have been done. 

After that, we checked also the stability of the results by performing further 
simulations with different values for the parameters of the Hamiltonian. We fixed the 
parameter D at the new value D = D* + 6D, where 5D = 0.008 is the confidence of 
the improved value D* = 0.641. Then, after having tuned the correct critical inverse 
temperature (3c [271 EHl ES], we performed Monte Carlo simulations on a lattice with 
length L = 96. We have done an average over 2-10^ Markov chains; for each chain 
we set different high-temperature initial conditions and, after the quench, we have done 
700 lattice sweeps using Metropolis algorithm. Performing the same data analysis that 
we have done in the previous sections we obtained an estimate of the dynamical critical 
exponent z = 2.022(6). Therefore, we argued that the uncertainty on the estimate of 
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Figure 4. Coefficient Cn for the fit Inx(i) = Inxo + plnt + Cut""! + Cait"''^ The 
dotted hne corresponds to the zero value. 



the critical value due to the numerical uncertainty of D* is 

5zc^\z{D*)- z{b)\^Q.m2{l). (22) 

Although this result is compatible with zero, we have to take into account its effects 
on the determination of the global incertitude of the value of the dynamical critical 
exponent z. Thus, we definitively predict 

2 = 2.020(6) [2], (23) 

where the error in brackets gives the variation of the estimate as D* varies within one 
error bar. 



3.2. Correlation Length 

Finally, using the correlation length, we have checked the results achieved in the last 
Section. In Figjsjwe plot ^{t, L) and ^{t,L)/L. Confronting this plot with Figjl] it is 
easy to see that when C,{t, L)/L < 0.22 there are not finite-size effects in the behaviour of 
susceptibihty x{t,L). In the subview of Figjsjthe gray stripe represents the data range 
into which susceptibility is starting to leave the infinite-volume curve. In the main 
graph we show the correlation length in log — log scale; the red straight line represents 
the infinite- volume behaviour according to the dynamical critical exponent that we have 
found in the last Section. Indeed, according to dynamical TZQ, in the thermodynamic 
limit, correlation length behaves as 

m = i.t'l^F^ (24) 

where is an analytic function that depends on the next-to-leading scaling fields. 

As shown by FigjS} also in this case the curves for different lattice size are parallels 
to the infinite- volume curve up to a value tmax{L) which increases as L^. 
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Figure 5. The behaviour of the correlation length: in the main plot ^(i, L) versus t 
in log-log scale; in the subview the rescaled ^(t, L)/L. 

4. Conclusions 

In this paper we liave studied tlie off-equilibrium purely relaxational dynamics in 
3d-Blume-Capel model. According to standard TZQ arguments applied to dynamics, 
the dynamical critical behaviour in such a system should belong to the same model- 
A dynamical universality class. If this description is correct, the dynamical critical 
exponent z is the same as in 3d-Ising model and the leading scaling corrections are 
controlled by the same TZQ operators that appear in the statics and therefore are 
characterized by the static correction-to-scaling exponents Ui = 0.84(4) and U2 = 
1.67(11). For the same reasons, in the case of improved Hamiltonians, leading scaling 
corrections should be also absent in dynamical quantities. Improved models are expected 
to provide the most precise estimates of universal dynamical quantities. For instance, 
in FSS studies at the critical point, corrections to scaling decay as L~^'^'^ in improved 
models, while in generic systems the approach to the thermodynamics limit is much 
slower, with corrections decaying as 

The main results of our analysis can be summarized as follows. 

We investigate the off-equilibrium relaxational dynamics in the 3d-Blume-Capel 
model. We started from disordered T = oo configurations and observed the relaxation 
at T = Tc. The results show that previous FT estimates of dynamical critical exponent 
z are perfectly consistent with our MC simulations that give the value z = 2.020(8). In 
the analysis particular attention has been taken to avoid finite-size corrections. 

Since the model is improved, we do not observe corrections proportional to t^"^^/^; 
instead our data show corrections that are proportional to t^^^^^. Moreover, we have 
been careful at new scaling corrections introduced by a new TZQ operator that appear 
in the dynamics and we have shown that it is possible to forget such a corrections 
throughout the data analysis. Indeed, the correlations of the order parameter at initial 
time surface (tuned by Tq"^) give rise to dynamical scaling corrections proportional to 
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t~^] nevertheless we can absorb such a corrections into t~^'^/^ . 

In conclusion, with this work, we have found an estimate of the dynamical 
critical exponent z for the three-dimensional Ising universality class that is perfectly 
in agreement with the FT predictions. 
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Appendix A. Metropolis Algorithm 

We have implemented the standard local Metropolis algorithm with the acceptance rate 

A = min[l,exp --H])] (A.l) 

where l-C and l-L correspond to the Hamiltonian evaluated for the proposal and for the 
given spin configuration, rispectively. The proposal is generated by choosing with the 
same probability one of the two new possible values of the spin at a single site x of the 
lattice. For example, if the initial value of the spin at site x is +1 then we randomly 
choose a new value from the set {0,-1}. Hence l-L' — l-L depends only on the value of the 
spin at the site x and at its neighbours y. We perform the single-site update sequentially, 
moving from one site to next-nearest-neighbour site in a checkerboard fashion. Defining 
the parity of a cubic lattice site as P{i,j, k) = (— 1)*+-'"'''^, we perform a sequential update 
first over all even sites then over all odd ones. The time unit is defined as a Monte Carlo 
sweep over all sites of the lattice. Moreover, the pseudorandom number generator is 
been implemented by using the 64-bit Mersenne Twister algorithm [5H 153] . 
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